This script:
meanbeta) in 8 focal regions (left and right SMG, left and
right STS, left and right V1, bilateral APC and RFC).It also contains the code needed to reproduce Figure 3 in the main
text of the paper. All instances of write_rds,
write_csv, etc have been commented out to avoid writing
over existing files. Uncomment those statements in order to reproduce
these steps and save (new) outputs.
univariate_data <- readRDS(here("outputs/univariate_data/study1_univariate_data_singlebetas.Rds"))
univariate_data$event <- relevel(as.factor(univariate_data$event), ref = "unexp")
univariate_summary_domain <- readRDS(here("outputs",
"univariate_descriptive_summaries",
"study1_univariate_summary_domain_allvoxelN.Rds"))
univariate_summary_task <- readRDS(here(
"outputs",
"univariate_descriptive_summaries",
"study1_univariate_summary_task_allvoxelN.Rds"))
focal_data <- univariate_data %>%
filter(focal_region == 1)
focal_summary_domain <- univariate_summary_domain %>%
filter(focal_region == 1)
focal_summary_task <- univariate_summary_task %>%
filter(focal_region == 1)
domain_specific_regions <- c("physics", "psychology")
domain_general_regions <- c("early_visual", "MD")
focal_data_100 <- focal_data %>%
filter(top_n_voxels == "100")
focal_summary_domain_100 <- focal_summary_domain %>%
filter(top_n_voxels == "100")
focal_summary_task_100 <- focal_summary_task %>%
filter(top_n_voxels == "100")
model_habituation <- lmer(data = focal_data_100 %>%
filter(ROI_category != "early_visual",
event != "fam",
ROI_category != "MD"), # cannot model MD regions because of non-independence of ROI selection (based on runs 2-4)
formula = meanbeta ~ extracted_run_number * event + (1|subjectID))
check_model(model_habituation)
plot(allEffects(model_habituation))
summary(model_habituation)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ extracted_run_number * event + (1 | subjectID)
## Data: focal_data_100 %>% filter(ROI_category != "early_visual", event !=
## "fam", ROI_category != "MD")
##
## REML criterion at convergence: 17809
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8142 -0.5931 -0.0313 0.5983 5.9668
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 0.2364 0.4862
## Residual 3.4448 1.8560
## Number of obs: 4352, groups: subjectID, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.181e-01 1.212e-01 1.600e+01 2.624 0.01843
## extracted_run_number1 5.679e-02 4.873e-02 4.328e+03 1.165 0.24396
## extracted_run_number2 -9.594e-02 4.873e-02 4.328e+03 -1.969 0.04904
## extracted_run_number3 7.969e-03 4.873e-02 4.328e+03 0.164 0.87010
## event1 7.583e-02 2.814e-02 4.328e+03 2.695 0.00706
## extracted_run_number1:event1 1.519e-01 4.873e-02 4.328e+03 3.117 0.00184
## extracted_run_number2:event1 -1.063e-02 4.873e-02 4.328e+03 -0.218 0.82738
## extracted_run_number3:event1 -6.623e-02 4.873e-02 4.328e+03 -1.359 0.17417
##
## (Intercept) *
## extracted_run_number1
## extracted_run_number2 *
## extracted_run_number3
## event1 **
## extracted_run_number1:event1 **
## extracted_run_number2:event1
## extracted_run_number3:event1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ext__1 ext__2 ext__3 event1 e__1:1 e__2:1
## extrctd_r_1 0.000
## extrctd_r_2 0.000 -0.333
## extrctd_r_3 0.000 -0.333 -0.333
## event1 0.000 0.000 0.000 0.000
## extrct__1:1 0.000 0.000 0.000 0.000 0.000
## extrct__2:1 0.000 0.000 0.000 0.000 0.000 -0.333
## extrct__3:1 0.000 0.000 0.000 0.000 0.000 -0.333 -0.333
model_habituation_results <- cbind(gen.m(model_habituation), gen.ci(model_habituation)[3:10,])
## Computing profile confidence intervals ...
model_habituation_results_byrun <- lsmeans(model_habituation, pairwise ~ event | extracted_run_number, pbkrtest.limit = 4352)$contrasts %>% as.data.frame()
DT::datatable(model_habituation_results, options = list(scrollX = TRUE))
DT::datatable(model_habituation_results_byrun, options = list(scrollX = TRUE))
# write_rds(model_habituation_results, path = here("outputs/univariate_results/habituation_over_runs_modelsummary.Rds"))
# write_rds(model_habituation_results_byrun, path = here("outputs/univariate_results/habituation_over_runs_perrun.Rds"))
regions <- levels(as.factor(focal_data$ROI_name_final))
plot_by_domain_run1 <- function(region, ymin, ymax) {
plotobject <- ggplot(data = focal_summary_domain_100 %>% filter(str_detect(ROI_name_final, region), extracted_run_number == "run1", event != "fam"),
aes(x = event, y = meanbeta, fill = domain)) +
geom_bar(stat = "identity", aes(alpha = event), colour = "black") +
geom_errorbar(
aes(ymin = meanbeta - se, ymax = meanbeta + se),
position = position_dodge(width = .9),
width = .2,
colour = "black"
) +
theme_cowplot(10) +
facet_wrap( ~ ROI_name_final + domain, nrow = 1) +
scale_fill_manual(values = c("#00AFBB", "#FC4E07")) +
ylab("Amplitude") +
xlab("Event") +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
),
legend.position = "none") +
coord_cartesian(ylim = c(ymin, ymax)) +
ggtitle(paste0("ROI:", region))
plotobject
}
This plot is a subset of Figure 3 in the paper.
SMG <- plot_by_domain_run1("supramarginal_", -1, 4)
STS <- plot_by_domain_run1("superiortemporal_", 0, 4)
early_vis <- plot_by_domain_run1("V1_", 0, 11)
MD_APC <- plot_by_domain_run1("antParietal_bilateral", -1, 1)
MD_RFC <- plot_by_domain_run1("precentral_inferiorfrontal_R", -1, 1)
(SMG + STS) / (early_vis + MD_APC + MD_RFC) + plot_layout(widths = c(1, 1, 2, 1))
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
regions <- levels(as.factor(focal_data$ROI_name_final))
plot_by_domain_all_runs <- function(region, ymin, ymax) {
plotobject <- ggplot(data = focal_summary_domain_100 %>% filter(str_detect(ROI_name_final, region)),
aes(x = event, y = meanbeta, fill = domain)) +
geom_bar(stat = "identity", aes(alpha = event), colour = "black") +
geom_errorbar(
aes(ymin = meanbeta - se, ymax = meanbeta + se),
position = position_dodge(width = .9),
width = .2,
colour = "black"
) +
theme_cowplot(10) +
facet_wrap( ~ ROI_name_final + top_n_voxels + domain + extracted_run_number, nrow = 1) +
scale_fill_manual(values = c("#00AFBB", "#FC4E07")) +
ylab("Average beta (amplitude)") +
xlab("Event") +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
),
legend.position = "none") +
coord_cartesian(ylim = c(ymin, ymax)) +
ggtitle(paste0("ROI:", region))
plotobject
}
plot_by_domain_all_runs("superiortemporal_L", -1, 3) | plot_by_domain_all_runs("superiortemporal_R", -1, 3)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_by_domain_all_runs("supramarginal_L", -1, 1) | plot_by_domain_all_runs("supramarginal_R", -1, 1)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_by_domain_all_runs("V1_L", -1, 7) | plot_by_domain_all_runs("V1_R", -1, 7)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_by_domain_run1("antParietal_bilateral", -1, 1) | plot_by_domain_run1("precentral_inferiorfrontal_R", -1, 1)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
regions <- levels(as.factor(focal_summary_task_100$ROI_name_final))
plot_univar_bytask <- function(region, ymin, ymax) {
plotobject <- ggplot(data = focal_summary_task_100 %>%
filter(str_detect(ROI_name_final, region),
extracted_run_number == "run1"),
aes(x = event, y = meanbeta, fill = task)) +
geom_bar(stat = "identity", aes(alpha = event), colour = "black") +
geom_errorbar(
aes(ymin = meanbeta - se, ymax = meanbeta + se),
position = position_dodge(width = .9),
width = .2,
colour = "black"
) +
# geom_point(data = focal_data_100 %>% filter(str_detect(ROI_name_final, region),
# extracted_run_number == "run1"),
# alpha = .1) +
# geom_line(data = focal_data_100 %>% filter(str_detect(ROI_name_final, region)),
# aes(group =
# subjectID),
# alpha = .1) +
theme_cowplot(10) +
facet_wrap( ~ ROI_name_final + factor(task, c("solidity", "support", "goal", "efficiency")), nrow = 1) +
scale_fill_manual(values = c("#deaf1f", "#f34b00", "#00aa8b", "#4f8e00")) +
ylab("Average beta (amplitude)") +
xlab("Event") +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
),
legend.position = "none") +
coord_cartesian(ylim = c(ymin, ymax)) +
ggtitle(paste0("ROI:", region))
plotobject
}
plot_univar_bytask("superiortemporal_L", -1, 3) | plot_univar_bytask("superiortemporal_R", -1, 3)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_univar_bytask("supramarginal_L", -1, 3) | plot_univar_bytask("supramarginal_R", -1, 3)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_univar_bytask("V1_L", -1, 7) | plot_univar_bytask("V1_R", -1, 7)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_univar_bytask("antParietal_bilateral", -1, 3) | plot_univar_bytask("precentral_inferiorfrontal_R", -1, 3)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
NOTE: Here we are using summed contrasts, which means that main effects are interpreted with respect to the grand mean.
effect_sizes_intercept_row <- data.frame(d = NA)
BF_intercept_row <- data.frame(BF = NA, BF_interpretation = NA)
# rownames(effect_sizes_intercept_row) <- "(Intercept)"
rownames(BF_intercept_row) <- "(Intercept)"
rownames_no_intercept <- c("event1", "domain1", "rep.L", "event:domain")
data_formodeling <- univariate_data %>%
filter(event != "fam")
data_formodeling %>%
group_by(top_n_voxels) %>%
summarise(n = n())
ROIs <- levels(as.factor(data_formodeling$ROI_name_final))
top_n_voxels <- unique(data_formodeling$top_n_voxels)
runs <- c("run1", "run2", "run3", "run4", "all_runs")
modelsummaries_init <- data.frame()
for (ROI in ROIs) {
curROI <- ROI
for (topN in top_n_voxels) {
cur_top_voxels_N <- topN
for (run in runs) {
curRun <- run
if (curRun != "all_runs") {
ROIdata <- data_formodeling %>%
filter(ROI_name_final == curROI,
top_n_voxels == topN,
extracted_run_number == curRun)
model <- lmer(data = ROIdata,
formula = meanbeta ~ event * domain + rep + (1|subjectID))
} else if (curRun == "all_runs") {
ROIdata <- data_formodeling %>%
filter(ROI_name_final == curROI,
top_n_voxels == topN)
model <- lmer(data = ROIdata,
formula = meanbeta ~ event * domain + rep + (1|subjectID))
}
model_just_maineffects <-
update(model, formula = ~ . - event:domain) # Without interaction term
model_just_domain <-
update(model_just_maineffects, formula = ~ . - event)
model_just_event <-
update(model_just_maineffects, formula = ~ . - domain)
model_no_rep <- update(model, formula = ~ . - rep)
BF_BIC_interaction <-
data.frame(BF = exp((
BIC(model_just_maineffects) - BIC(model)
) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_just_maineffects) - BIC(model)
) / 2))) # BICs to Bayes factor
BF_BIC_event <-
data.frame(BF = exp((
BIC(model_just_domain) - BIC(model_just_maineffects)
) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_just_domain) - BIC(model_just_maineffects)
) / 2)))
BF_BIC_domain <-
data.frame(BF = exp((
BIC(model_just_event) - BIC(model_just_maineffects)
) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_just_event) - BIC(model_just_maineffects)
) / 2)))
BF_BIC_rep <-
data.frame(BF = exp((BIC(model_no_rep) - BIC(model)) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_no_rep) - BIC(model)
) / 2)))
BF_initial <- rbind(BF_BIC_event,
BF_BIC_domain,
BF_BIC_rep,
BF_BIC_interaction)
rownames(BF_initial) <- rownames_no_intercept
BFs <- rbind(BF_intercept_row,
BF_initial)
effect_sizes <- standardize_parameters(model) %>%
as.data.frame() %>%
select(Std_Coefficient)
if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
ROI_category <- "domain_specific"
} else {
ROI_category <- "domain_general"
}
ROI_focal <- ROIdata$focal_region[1]
modelsummary <-
cbind(
cur_top_voxels_N,
curRun,
ROI_category,
ROI_focal,
curROI,
rownames(summary(model)$coefficients),
summary(model)$coefficients,
confint.merMod(model)[3:7, ],
effect_sizes,
BFs
) %>%
as.data.frame()
colnames(modelsummary) <-
c(
"top_n_voxels",
"extracted_run_number",
"domain",
"focal_region",
"ROI",
"effect",
"B",
"SE",
"df",
"t",
"p",
"CI_95_lower",
"CI_95_upper",
"B_standardized",
"BF",
"BF_interpretation"
)
modelsummaries_init <-
rbind(modelsummaries_init, modelsummary)
}
}
}
modelsummaries <- modelsummaries_init %>%
mutate_at(c(7:15), as.numeric) %>%
mutate(star = as.factor(
case_when(p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .1 ~ "~",
TRUE ~ " ")
)) %>%
mutate(p = round(p, digits = 3))
modelsummaries$effect <-
as.factor(modelsummaries$effect)
levels(modelsummaries$effect) <-
c("intercept", "domain", "event", "event:domain", "rep")
rownames(modelsummaries) <- NULL
# write_rds(modelsummaries, here("outputs", "univariate_results", "modelsummaries_allregions_alltopN.Rds"))
focal_modelsummaries <- read_rds(here("outputs/univariate_results/modelsummaries_allregions_alltopN.Rds")) %>%
filter(focal_region == 1)
DT::datatable(focal_modelsummaries %>% filter(top_n_voxels == "100", focal_region == 1) %>% arrange(effect), options = list(scrollX = TRUE,
pageLength = 100))
# filter(extracted_run_number == "all_runs" | extracted_run_number == "run1") %>%
LSMG_model <- lmer(data = univariate_data %>%
filter(top_n_voxels == "100",
event != "fam",
extracted_run_number == "run1",
ROI_name_final == "supramarginal_L"),
formula = meanbeta ~ event * domain + rep + (1|subjectID))
summary(LSMG_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + rep + (1 | subjectID)
## Data: univariate_data %>% filter(top_n_voxels == "100", event != "fam",
## extracted_run_number == "run1", ROI_name_final == "supramarginal_L")
##
## REML criterion at convergence: 928.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.85161 -0.63506 -0.00846 0.60802 2.64499
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 0.5561 0.7457
## Residual 1.5254 1.2351
## Number of obs: 272, groups: subjectID, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.01515 0.19576 16.00000 0.077 0.939275
## event1 0.26528 0.07489 251.00000 3.542 0.000473 ***
## domain1 0.38079 0.07489 251.00000 5.085 7.2e-07 ***
## rep.L -0.07030 0.10591 251.00000 -0.664 0.507425
## event1:domain1 0.25040 0.07489 251.00000 3.344 0.000953 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 rep.L
## event1 0.000
## domain1 0.000 0.000
## rep.L 0.000 0.000 0.000
## event1:dmn1 0.000 0.000 0.000 0.000
LSMG_results_bydomain_top100 <- lsmeans(LSMG_model, pairwise ~ event | domain)$contrasts %>% as.data.frame()
DT::datatable(LSMG_results_bydomain_top100, options = list(scrollX = TRUE))
# write_rds(LSMG_results_bydomain_top100, path = here("outputs/univariate_results/LSMG_results_bydomain_top100.Rds"))
RSMG_model <- lmer(data = univariate_data %>%
filter(top_n_voxels == "100",
event != "fam",
extracted_run_number == "run1",
ROI_name_final == "supramarginal_R"),
formula = meanbeta ~ event * domain + rep + (1|subjectID))
summary(RSMG_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + rep + (1 | subjectID)
## Data: univariate_data %>% filter(top_n_voxels == "100", event != "fam",
## extracted_run_number == "run1", ROI_name_final == "supramarginal_R")
##
## REML criterion at convergence: 1047.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2569 -0.5727 0.0203 0.6278 2.5903
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 1.517 1.232
## Residual 2.312 1.520
## Number of obs: 272, groups: subjectID, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.01345 0.31266 16.00000 -0.043 0.96621
## event1 0.21889 0.09219 251.00000 2.374 0.01834 *
## domain1 0.28691 0.09219 251.00000 3.112 0.00207 **
## rep.L -0.36357 0.13038 251.00000 -2.789 0.00570 **
## event1:domain1 0.05727 0.09219 251.00000 0.621 0.53505
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 rep.L
## event1 0.000
## domain1 0.000 0.000
## rep.L 0.000 0.000 0.000
## event1:dmn1 0.000 0.000 0.000 0.000
lsmeans(RSMG_model, pairwise ~ event | domain)$contrasts
## domain = physics:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.552 0.261 251 2.118 0.0352
##
## domain = psychology:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.323 0.261 251 1.240 0.2163
##
## Results are averaged over the levels of: rep
## Degrees-of-freedom method: kenward-roger
LSTS_model <- lmer(data = univariate_data %>%
filter(top_n_voxels == "100",
event != "fam",
extracted_run_number == "run1",
ROI_name_final == "superiortemporal_L"),
formula = meanbeta ~ event * domain + rep + (1|subjectID))
summary(LSTS_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + rep + (1 | subjectID)
## Data: univariate_data %>% filter(top_n_voxels == "100", event != "fam",
## extracted_run_number == "run1", ROI_name_final == "superiortemporal_L")
##
## REML criterion at convergence: 1051.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4985 -0.5627 -0.0406 0.6461 2.9036
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 1.773 1.331
## Residual 2.327 1.525
## Number of obs: 272, groups: subjectID, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.54190 0.33591 16.00000 1.613 0.126244
## event1 0.16841 0.09249 251.00000 1.821 0.069839 .
## domain1 -0.32878 0.09249 251.00000 -3.555 0.000452 ***
## rep.L 0.13244 0.13081 251.00000 1.012 0.312286
## event1:domain1 0.04951 0.09249 251.00000 0.535 0.592897
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 rep.L
## event1 0.000
## domain1 0.000 0.000
## rep.L 0.000 0.000 0.000
## event1:dmn1 0.000 0.000 0.000 0.000
lsmeans(LSTS_model, pairwise ~ event | domain)$contrasts
## domain = physics:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.436 0.262 251 1.666 0.0970
##
## domain = psychology:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.238 0.262 251 0.909 0.3643
##
## Results are averaged over the levels of: rep
## Degrees-of-freedom method: kenward-roger
RSTS_model <- lmer(data = univariate_data %>%
filter(top_n_voxels == "100",
event != "fam",
extracted_run_number == "run1",
ROI_name_final == "superiortemporal_R"),
formula = meanbeta ~ event * domain + rep + (1|subjectID))
summary(RSTS_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + rep + (1 | subjectID)
## Data: univariate_data %>% filter(top_n_voxels == "100", event != "fam",
## extracted_run_number == "run1", ROI_name_final == "superiortemporal_R")
##
## REML criterion at convergence: 1015.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.16150 -0.59428 -0.03496 0.59256 2.87577
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 2.143 1.464
## Residual 1.992 1.411
## Number of obs: 272, groups: subjectID, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.95600 0.36522 16.00000 2.618 0.018662 *
## event1 0.25840 0.08557 251.00000 3.020 0.002792 **
## domain1 -0.30885 0.08557 251.00000 -3.609 0.000371 ***
## rep.L -0.22514 0.12102 251.00000 -1.860 0.064000 .
## event1:domain1 -0.05083 0.08557 251.00000 -0.594 0.553030
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 rep.L
## event1 0.000
## domain1 0.000 0.000
## rep.L 0.000 0.000 0.000
## event1:dmn1 0.000 0.000 0.000 0.000
lsmeans(RSTS_model, pairwise ~ event | domain)$contrasts
## domain = physics:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.415 0.242 251 1.715 0.0876
##
## domain = psychology:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.618 0.242 251 2.555 0.0112
##
## Results are averaged over the levels of: rep
## Degrees-of-freedom method: kenward-roger
NOTE: Here we are using summed contrasts, which means that main effects are interpreted with respect to the grand mean.
options(contrasts = c("contr.sum", "contr.poly"))
effect_sizes_intercept_row <- data.frame(d = NA)
BF_intercept_row <- data.frame(BF = NA, BF_interpretation = NA)
# rownames(effect_sizes_intercept_row) <- "(Intercept)"
rownames(BF_intercept_row) <- "(Intercept)"
rownames_no_intercept <- c("event1", "rep.L")
data_formodeling <- univariate_data %>%
filter(event != "fam",
top_n_voxels == 100)
data_formodeling %>%
group_by(top_n_voxels) %>%
summarise(n = n())
ROIs <- levels(as.factor(data_formodeling$ROI_name_final))
top_n_voxels <- unique(data_formodeling$top_n_voxels)
runs <- c("run1", "all_runs")
tasks <- unique(data_formodeling$task)
ROIs <- unique(data_formodeling$ROI_name_final)
data_formodeling$event <- relevel(data_formodeling$event, ref = "unexp")
modelsummaries_pertask_init <- data.frame()
for (ROI in ROIs) {
curROI <- ROI
for (task in tasks) {
cur_task <- task
for (run in runs) {
if (run == "run1") {
ROIdata <-
data_formodeling %>% filter(ROI_name_final == curROI,
task == cur_task,
extracted_run_number == "run1")
} else if (run == "all_runs")
ROIdata <-
data_formodeling %>% filter(ROI_name_final == curROI,
task == cur_task)
model <- lmer(data = ROIdata,
formula = meanbeta ~ event + rep + (1 | subjectID)) # had to remove run random int (convergence))
focal_region <- ROIdata$focal_region[1]
task_domain <- as.character(ROIdata$domain[1])
ROI_category <- ROIdata$ROI_category[1]
model_null <-
update(model, formula = ~ . - event) # Without interaction term
model_norep <-
update(model, formula = ~ . - rep) # Without interaction term
BF_BIC_event <-
data.frame(BF = exp((BIC(model_null) - BIC(model)) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_null) - BIC(model)
) / 2))) # BICs to Bayes factor
BF_BIC_rep <-
data.frame(BF = exp((BIC(model_norep) - BIC(model)) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_norep) - BIC(model)
) / 2))) # BICs to Bayes factor
# BF_initial <- rbind(BF_BIC_event)
BFs <- rbind(BF_intercept_row,
BF_BIC_event,
BF_BIC_rep)
rownames(BFs)[2:3] <- rownames_no_intercept
effect_sizes <- standardize_parameters(model) %>%
as.data.frame() %>%
select(Std_Coefficient)
cur_top_voxels_N <- 100
modelsummary <-
cbind(
cur_top_voxels_N,
run,
ROI_category,
focal_region,
curROI,
cur_task,
rownames(summary(model)$coefficients),
summary(model)$coefficients,
confint.merMod(model)[3:5, ],
effect_sizes,
BFs
) %>%
as.data.frame()
colnames(modelsummary) <-
c(
"top_n_voxels",
"extracted_run_number",
"domain",
"focal_region",
"ROI",
"task",
"effect",
"B",
"SE",
"df",
"t",
"p",
"CI_95_lower",
"CI_95_upper",
"B_standardized",
"BF",
"BF_interpretation"
)
modelsummaries_pertask_init <-
rbind(modelsummaries_pertask_init, modelsummary)
}
}
}
modelsummaries_pertask <- modelsummaries_pertask_init %>%
mutate_at(c(8:16), as.numeric) %>%
mutate(star = as.factor(
case_when(p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .1 ~ "~",
TRUE ~ " ")
)) %>%
mutate(p = round(p, digits = 3))
modelsummaries_pertask$effect <-
as.factor(modelsummaries_pertask$effect)
levels(modelsummaries_pertask$effect) <-
c("intercept", "event", "rep")
rownames(modelsummaries_pertask) <- NULL
# write_rds(modelsummaries_pertask, here("outputs", "univariate_results", "modelsummaries_allregions_pertask.Rds"))
modelsummaries_pertask <- read_rds(here("outputs", "univariate_results", "modelsummaries_allregions_pertask.Rds"))
focal_modelsummaries_pertask_forplotting <- modelsummaries_pertask %>%
filter(top_n_voxels == "100",
effect == "event",
extracted_run_number == "run1",
focal_region == 1)
focal_modelsummaries_pertask_forplotting$ROI <- factor(focal_modelsummaries_pertask_forplotting$ROI, levels = c("supramarginal_L", "supramarginal_R", "superiortemporal_L", "superiortemporal_R", "antParietal_bilateral", "precentral_inferiorfrontal_R", "V1_L", "V1_R"))
pertask_b <- ggplot(focal_modelsummaries_pertask_forplotting,
aes(ROI, B, colour = task, group = task)) +
geom_point() +
geom_errorbar(
aes(ymin = B - SE, ymax = B + SE,
colour = task),
position = position_dodge(width = .9),
width = .2
) +
# geom_line() +
geom_hline(yintercept = 0) +
facet_wrap( ~ task, nrow = 1) +
scale_colour_manual(values = c("#deaf1f", "#f34b00", "#00aa8b", "#4f8e00")) +
ylab("B and SE") +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
))
coord_cartesian(ylim = c(-.5, .7))
## <ggproto object: Class CoordCartesian, Coord, gg>
## aspect: function
## backtransform_range: function
## clip: on
## default: FALSE
## distance: function
## expand: TRUE
## is_free: function
## is_linear: function
## labels: function
## limits: list
## modify_scales: function
## range: function
## render_axis_h: function
## render_axis_v: function
## render_bg: function
## render_fg: function
## setup_data: function
## setup_layout: function
## setup_panel_guides: function
## setup_panel_params: function
## setup_params: function
## train_panel_guides: function
## transform: function
## super: <ggproto object: Class CoordCartesian, Coord, gg>
perregion_b <- ggplot(focal_modelsummaries_pertask_forplotting,
aes(task, B, colour = task, group = task)) +
geom_point() +
geom_errorbar(
aes(ymin = B - SE, ymax = B + SE,
colour = task),
position = position_dodge(width = .9),
width = .2
) +
# geom_line() +
geom_hline(yintercept = 0) +
facet_wrap( ~ ROI, nrow = 1) +
scale_colour_manual(values = c("#deaf1f", "#f34b00", "#00aa8b", "#4f8e00")) +
ylab("B and SE") +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
)) +
coord_cartesian(ylim = c(-.5, .7))
pertask_b
perregion_b
focal_modelsummaries$top_n_voxels <- factor(focal_modelsummaries$top_n_voxels, levels = c("10", "50", "100", "150", "200", "250", "300"))
focal_modelsummaries$extracted_run_number <- factor(focal_modelsummaries$extracted_run_number, levels = c("run1", "run2", "run3", "run4", "all_runs"))
ggplot(focal_modelsummaries %>%
filter(extracted_run_number != "all_runs",
ROI != "antParietal_bilateral",
ROI != "precentral_inferiorfrontal_R",
effect != "intercept",
effect != "rep"),
aes(extracted_run_number, B_standardized, colour = top_n_voxels)) +
geom_point() +
geom_line(aes(group = top_n_voxels)) +
ylab("Standardized Beta") +
facet_wrap(~ effect + factor(ROI, levels = c("supramarginal_L", "supramarginal_R", "superiortemporal_L", "superiortemporal_R", "V1_L", "V1_R")), nrow = 3) +
# coord_cartesian(ylim = c(-.75, 0.75)) +
geom_hline(yintercept = 0) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
nonfocal_modelsummaries <- read_rds(here("outputs/univariate_results/modelsummaries_allregions_alltopN.Rds")) %>%
filter(focal_region == 0)
DT::datatable(nonfocal_modelsummaries %>%
filter(top_n_voxels == "100") %>%
arrange(effect),
options = list(scrollX = TRUE,
pageLength = 100))
Here we repeat the MVPA manyregions analysis, whose results are difficult to interpret given the lack of reliable multivariate signal towards unexpected events, despite clear univariate effects. NOTE: Here we are using treatment contrasts, which means that main effects are interpreted with respect to the reference level.
options(contrasts = c("contr.treatment", "contr.poly"))
data_for_modeling <- univariate_data %>%
filter(event != "fam",
top_n_voxels == "100") %>%
filter(extracted_run_number == "run1") %>%
filter(!str_detect(ROI_name_final, "bilateral"))
data_for_modeling$event <- relevel(data_for_modeling$event, ref = "exp")
data_for_modeling$domain <- relevel(data_for_modeling$domain, ref = "psychology")
ROIs <- levels(as.factor(data_for_modeling$ROI_name_final))
# effect_sizes_intercept_row <- data.frame(d = NA)
rownames_no_intercept <- c("event1")
modelsummaries_init <- data.frame()
for (ROI in ROIs) {
curROI <- ROI
# FIRST COMPUTE EVENT EFFECT SIZES, FOR EACH DOMAIN
physics_data <- data_for_modeling %>%
filter(ROI_name_final == curROI,
domain == "physics")
psychology_data <- data_for_modeling %>%
filter(ROI_name_final == curROI,
domain == "psychology")
# had to remove random intercept for runs (convergence))
physics_model <- lmer(data = physics_data,
formula = meanbeta ~ event + (1|subjectID))
psychology_model <- lmer(data = psychology_data,
formula = meanbeta ~ event + (1|subjectID))
physics_event_effect_size <- standardize_parameters(physics_model) %>%
as.data.frame() %>%
select(Std_Coefficient) %>%
slice(2)
psychology_event_effect_size <- standardize_parameters(psychology_model) %>%
as.data.frame() %>%
select(Std_Coefficient) %>%
slice(2)
# SECOND COMPUTE DOMAIN EFFECT SIZES, FOR EACH EVENT
expected_data <- data_for_modeling %>%
filter(ROI_name_final == curROI,
event == "exp")
unexpected_data <- data_for_modeling %>%
filter(ROI_name_final == curROI,
event == "unexp")
# had to remove random intercept for runs (convergence))
exp_model <- lmer(data = expected_data,
formula = meanbeta ~ domain + (1 | subjectID))
unexp_model <- lmer(data = unexpected_data,
formula = meanbeta ~ domain + (1 | subjectID))
exp_domain_effect_size <-standardize_parameters(exp_model) %>%
as.data.frame() %>%
select(Std_Coefficient) %>%
slice(2)
unexp_domain_effect_size <- standardize_parameters(unexp_model) %>%
as.data.frame() %>%
select(Std_Coefficient) %>%
slice(2)
if (physics_data$ROI_category[1] == "MD" | physics_data$ROI_category[1] == "early_visual") {
ROI_domain <- "general"
} else {
ROI_domain <- "specific"
}
modelsummary <-
cbind(
ROI_domain,
physics_data$ROI_category[1],
curROI,
physics_event_effect_size,
psychology_event_effect_size,
exp_domain_effect_size,
unexp_domain_effect_size
) %>%
as.data.frame()
modelsummaries_init <-
rbind(modelsummaries_init, modelsummary)
}
colnames(modelsummaries_init) <-
c("ROI_domain",
"ROI_category",
"ROI",
"physics_event_effect_size",
"psychology_event_effect_size",
"exp_domain_effect_size",
"unexp_domain_effect_size"
)
rownames(modelsummaries_init) <- NULL
# write_rds(modelsummaries_init, here("outputs/univariate_results/univar_manyregions_results.Rds"))
univar_manyregions_results <- read_rds( here("outputs/univariate_results/univar_manyregions_results.Rds"))
plot1 <- ggplot(univar_manyregions_results, aes(psychology_event_effect_size, physics_event_effect_size)) +
geom_hline(yintercept = 0, linetype = 'dotted') +
geom_vline(xintercept = 0, linetype = 'dotted') +
coord_cartesian(xlim = c(-0.2, .75), ylim = c(-0.2, .75)) + geom_point(aes(colour = ROI_category), size = 3) +
geom_smooth(method = "lm") +
facet_wrap(~ROI_domain, scales = "free") +
# coord_cartesian(xlim = c(-1.5,1.5), ylim = c(-1.5,1.5)) +
scale_colour_manual(values = c("#fb00d3", "#f8bf00", "#4894ff", "#f71d00")) +
xlab("VOE Effect Size (Standardized Beta), Physics Events") +
ylab("VOE Effect Size (Standardized Beta), Psychology Events") +
ggtitle(label = "Exp 1 VOE effects for each domain across regions") +
geom_text_repel(aes(label = ROI), size = 3) +
theme(axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
legend.position = "none")
plot2 <- ggplot(univar_manyregions_results, aes(exp_domain_effect_size, unexp_domain_effect_size)) +
geom_hline(yintercept = 0, linetype = 'dotted') +
geom_vline(xintercept = 0, linetype = 'dotted') +
coord_cartesian(xlim = c(-1.2, 1.2), ylim = c(-1.2, 1.2)) +
geom_point(aes(colour = ROI_category), size = 3) +
geom_smooth(method = "lm") +
facet_wrap(~ROI_domain, scales = "free") +
# coord_cartesian(xlim = c(-1.5,1.5), ylim = c(-1.5,1.5)) +
scale_colour_manual(values = c("#fb00d3", "#f8bf00", "#4894ff", "#f71d00")) +
xlab("Domain Effect Size (Standardized Beta), Expected Events") +
ylab("Domain Effect Size (Standardized Beta), Unexpected Events") +
ggtitle(label = "Exp 1 Domain effects for each event across regions") +
geom_text_repel(aes(label = ROI), size = 3) +
theme(axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
legend.position = "none")
plot1
## `geom_smooth()` using formula = 'y ~ x'
plot2
## `geom_smooth()` using formula = 'y ~ x'
# instead of testing whether the linear relationship between x and y is 0,
# we test for independence instead. H0 is that x and y are independent; F_{XY}(x,y) = F_X(x) F_Y(y).
observed_cors_ind <- univar_manyregions_results %>%
group_by(ROI_domain) %>%
summarise(
cor_domain_within_event =
np.cor.test(
exp_domain_effect_size,
unexp_domain_effect_size,
alternative = "two.sided",
independent = TRUE
)$estimate,
p_domain_within_event =
np.cor.test(
exp_domain_effect_size,
unexp_domain_effect_size,
alternative = "two.sided",
independent = TRUE
)$p.value,
cor_event_within_domain =
np.cor.test(
physics_event_effect_size,
psychology_event_effect_size,
alternative = "two.sided",
independent = TRUE
)$estimate,
p_event_within_domain =
np.cor.test(
physics_event_effect_size,
psychology_event_effect_size,
alternative = "two.sided",
independent = TRUE
)$p.value
)
observed_cors_ind
# write_rds(observed_cors_ind, path = here("outputs/univariate_results/study1_univar_manyregions_observed_cor.Rds"))
set.seed(2020)
## First define a function to work out the difference of corrs:
diff_corr <- function(data, indices) {
data <- data[indices,]
cor1 <-
np.cor.test(DS_results$exp_domain_effect_size,
DS_results$unexp_domain_effect_size,
alternative = "two.sided",
independent = TRUE,
parallel = TRUE)$estimate
cor2 <-
np.cor.test(
data$psychology_event_effect_size,
data$physics_event_effect_size,
alternative = "two.sided",
independent = TRUE,
parallel = TRUE)$estimate
return(cor1 - cor2)
}
DS_results <- univar_manyregions_results %>%
filter(ROI_domain == "specific")
DG_results <- univar_manyregions_results %>%
filter(ROI_domain == "general")
# doMC::registerDoMC(cores = 2) # for running in parallel
# Then apply a bootstrap procedure with 4000 draws (uncomment to reproduce):
# res_boot_DS <- boot(
# data = DS_results,
# R = 4000,
# statistic = diff_corr,
# stype = "i"
# )
# saveRDS(res_boot_DS, here("outputs/univariate_results/study1_univariate_dsregions_4000perms_confirmatory.Rds"))
res_boot_DS <- read_rds(here("outputs/univariate_results/study1_univariate_dsregions_4000perms_confirmatory.Rds"))
## Retrieve the empirical 95% confidence interval:
ds_results <- append(boot.ci(res_boot_DS, type = "perc", conf = 0.95),boot.pval(res_boot_DS))
# saveRDS(ds_results, here("outputs/univariate_results/study1_univariate_dsregions_summary.Rds"))
## Then apply a bootstrap procedure with 4000 draws (uncomment to reproduce):
# res_boot_DG <- boot(data = DG_results,
# R = 4000,
# statistic = diff_corr,
# stype = "i")
# saveRDS(res_boot_DG, here("outputs/univariate_results/study1_univariate_dgregions_4000perms_confirmatory.Rds"))
res_boot_DG <- readRDS(here("outputs/univariate_results/study1_univariate_dgregions_4000perms_confirmatory.Rds"))
plot(res_boot_DG)
## Retrieve the empirical 95% confidence interval:
dg_results <- append(boot.ci(res_boot_DG, type = "perc", conf = 0.95),boot.pval(res_boot_DG))
# saveRDS(dg_results, here("outputs/univariate_results/study1_univariate_dgregions_summary.Rds"))